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This paper proposes a new, potentially useful, way in which to characterise the de- 
gree of stochasticity exhibited by orbits in a fixed time-independent galactic potential. 
This approach differs from earlier work involving the computation of Liapounov ex- 
ponents in two ways, namely (1) by focusing on the statistical properties of ensembles 
of trajectories, rather than individual orbits, and (2) by restricting attention to the 
properties of these ensembles over time scales shorter than the age of the Universe, tn- 
For many potentials, generic ensembles of initial conditions corresponding to stochas- 
tic orbits will evolve relatively quickly towards a time-independent invariant measure 
r, which is arguably the natural unit to consider if one is interested in self-consistent 
equilibria. The basic idea proposed here is to compute short time Liapounov charac- 
teristic numbers xi^'t) over time intervals At for orbits in an ensemble that samples 
this invariant measure, and to analyse the overall distribution of these x's- This is 
done in detail for one model potential, namely the sixth order truncation of the Toda 
lattice potential. One especially significant conclusion is that time averages and en- 
semble averages coincide, so that the form of the distribution of short time x('^^)'s 
for such an ensemble is actually encoded in the calculation of x(^) for a single orbit 
over long times t 3> Ih- The distribution of short time x's is analysed as a function 
of the energy E of orbits in the ensemble and the length of the short time sampling 
interval At. For relatively high energies, the distribution is essentially Gaussian, the 
dispersion decreasing with time as t~^, with an exponent < p < 1/2 that depends 
on the energy E. 
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1. MOTIVATION 



Dating back to the pioneering work of Henon and Heiles (1964) thirty years ago, con- 
siderable attention has focused on the study of orbits in nonintegral potentials which 
are thought to reflect the bulk mass distribution of a galaxy (see, e.g.. Martinet & 
Pfenniger 1987 or Contopoulos & Grosbol 1989, and references cited therein). The 
Hamiltonian systems associated with these potentials differ from integrable Hamilto- 
nian systems in that they may admit both regular and stochastic orbits. 

In order to characterise the stochastic orbits, it is useful to introduce the concept 
of a Liapounov exponent. One way in which to define such Liapounov exponents (cf. 
Pesin 1977) is in a fashion related directly to the Kolmogorov entropy (cf. Chirikov 
1979). However, it quite natural physically and more convenient computationally 
to define these exponents instead in terms of the stability of individual stochastic 
trajectories (Bennetin et al 1976) with respect to small perturbations. Specifically, 
Liapounov exponents may be defined by the prescription 



where 5r{0) and 5r{t) denote respectively the configuration space deviations of two 
nearby orbits at times and t. Of particular interest is the maximal Liapounov 
exponent corresponding to the most unstable perturbation, which can be computed 
numerically by selecting 5r{0) at random. In everything that follows, the words 
"Liapounov exponent" will refer to this maximal Liapounov exponent. 

This prescription involves a formal i — > oo limit. Clearly, however, one cannot 
integrate forever, so that, as a practical matter, one is restricted to computing finite 
times estimates of X. The crucial question then is: For how long must one integrate 
to obtain reasonable estimates? 

The precise answer to this question depends on the particular form of the potential. 
However, one knows from experience (cf. Contopoulos and Barbanis 1989) that, as 
a general rule, one must typically integrate over extremely long times, say t > IOHd, 
where denotes a characteristic dynamical, or crossing, time. Unfortunately, though, 
t/) ~ 10^ yr for a galaxy like the Milky Way, so that this entails an integration for a 
period of time that is orders of magnitude longer than the age of the Universe. For 
this reason, the Liapounov exponent X docs not provide a useful characterisation of 
the instability of individual trajectories on short time scales < tn- 

Motivated by this realisation, a number of different workers, notably Udry and 
Pfenniger (1988), have sought instead to define shorter time analogues of the Lia- 
pounov exponent which could provide a characterisation of the stochasticity of indi- 
vidual orbits on time scales < tn- This is an extremely interesting idea. However, 
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these authors have not yet fully explored the potential implications of such an ap- 
proach. 

Most of the work hitherto on the problem of stochasticity in galactic potentials has 
been formulated in terms of the properties of individual orbits. However, even though 
it is convenient computationally to study such individual trajectories, there are solid 
physical motivations for working instead with ensembles of orbits. Perhaps the most 
obvious point is that individual orbits of stars are not accessible observationally: all 
that one can detect is an instantaneous snapshot of the overall surface brightness of 
a galaxy which, one hopes, traces the distribution of mass. Another obvious point is 
that one is really compelled to consider ensembles of orbits, rather than individual 
trajectories, if one is interested in incorporating stochastic orbits into a self-consistent 
model using Schwarzschild's (1979) method, or any variant thereof. 

The object of this paper is to propose an alternative approach to the analysis of 
stochastic orbits which confronts these basic difficulties. Specifically, the aim is to 
provide useful characterisations of stochasticity appropriate for ensembles of orbits, 
focusing exclusively on relatively short time scales < tn- 

A consideration of stochasticity for individual orbits on short time scales is a very 
complicated proposition and, for this reason, has been avoided by most nonlinear 
dynamicists. However, even though the transient behavior of individual orbits is quite 
complex, the statistical properties of ensembles of orbits can, nevertheless, exhibit 
striking regularities which facilitate a relatively simple characterisation. 

Three basic facts, described in detail below, facilitate this general approach. The 
first is that, for many potentials, including the model system analysed in this paper 
(Kandrup & Mahon 1994), generic ensembles of initial conditions corresponding to 
stochastic orbits will evolve on a relatively short time scale towards a particular dis- 
tribution which is time-independent or, if not strictly time- independent, only exhibits 
subsequent variability on very long time scales. The approach towards a strictly time- 
independent invariant measure is in fact guaranteed for Hamiltonian systems in which 
the stochastic orbits are ergodic, although there is no reason theoretically to expect 
that the time scale associated with this approach is short. 

The second fact is that, in many contexts, it is this invariant, or near- invariant, 
measure which constitutes the natural object in terms of which to analyse the statisti- 
cal properties of stochasticity. In particular, it is reasonable to consider an ensemble 
of stochastic orbits that samples the invariant measure, and then compute the Li- 
apounov characteristic numbers for the orbits in that ensemble over time intervals 
At < tn- The statistical properties of the resulting distribution of Liapounov char- 
acteristic numbers then provides a useful characterisation of the overall degree of 
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stochasticity for the ensemble. 

The final fact is that, for an ensemble that samples the invariant measure, ensemble 
averages and time averages may be equivalent, so that all the statistical information 
about stochasticity appropriate for such ensembles can actually be contained within 
the standard calculation of a Liapounov characteristic number X{t) for a single orbit 
over very long times t. For the case of a true invariant measure, this equivalence 
follows from the ergodic theorem. 

These three facts indicate that it is especially natural to compute short time 
Liapounov characteristic numbers for the particular case of ensembles of orbits that 
sample the invariant measure. However, one can equally well choose to compute 
these characteristic numbers for other ensembles. Thus, for example, one might focus 
upon an ensemble initially localised in some small phase space region which eventually 
evolves towards the invariant measure, and then compute the Liapounov characteristic 
numbers for the orbits in that ensemble. Such a computation is potentially interesting 
in that it permits a characterisation of the overall stochasticity of a "nonequilibrium" 
ensemble as it "evolves towards equilibrium," which can facilitate useful insights into 
the connection between stochasticity and relaxation. 

Section 2 of this paper outlines the general tact to be used in providing short time 
characterisations of stochasticity for ensembles of orbits. The remaining sections then 
focus on an implementation of this approach for one particular model, namely the 
sixth order truncation of the Toda (1967a,b) lattice potential. This specific form was 
not selected because of the expectation that it provides an especially good approxi- 
mation to the gravitational potential of any particular class of galaxies. Rather, in 
the spirit of the pioneering work of Henon and Heiles (1964), the aim has been to 
choose a potential which, in the sense described in Section 3, appears to be generic; 
and to use that potential as a testing ground for various ideas. 

Section 3 of the paper first demonstrates the concrete sense in which the informa- 
tion about short time measures of stochasticity is actually encoded in the long time 
calculations of X{t), and then analyses the form of the distribution of short time Lia- 
pounov characteristic numbers associated with an ensemble of orbits that samples the 
invariant measure. Section 4 next outlines one concrete prescription for generating 
ensembles that sample this invariant measure by evolving initially localised ensembles 
of initial conditions, and then analyses the distribution of Liapounov characteristic 
numbers associated with these initially localised ensembles. Section 5 concludes by 
examining the form of these distributions as a function of the sampling interval At. 
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2. TRANSIENT DYNAMICS OF AN ENSEMBLE OF STOCHASTIC 
ORBITS 



The standard formulation of nonlinear dynamics, as applied to galactic dynamics, en- 
tails a consideration of asymptotic orbital dynamics, focusing explicitly on the asymp- 
totic, long time behavior of individual trajectories. The object here is to propose a 
new approach, transient ensemble dynamics. In this alternative approach, the princi- 
pal focus is on the statistical properties of ensembles of orbits, rather than the details 
of individual trajectories. Moreover, the analysis is restricted entirely to properties 
of these orbits on short time scales < tn- In this approach, asymptotic calculations 
over extremely long time intervals are eschewed as unphysical, except to the extent 
that they can provide information about the behavior of ensembles of orbits on short, 
astrophysically relevant, time scales. 

As noted already, the transient behavior of individual orbits can be extremely 
difficult to characterise, so that a description of stochasticity in terms of transient or- 
bital dynamics would typically prove quite complicated. However, there are intuitive 
reasons to believe that, nevertheless, the physically relevant aspects of the short time 
behavior can be characterised relatively easily if one chooses instead to focus on the 
statistical properties of ensembles of orbits. 

The key feature entering into the proposed statistical description of stochastic 
orbits is the notion of an invariant measure. Mathematically, the invariant mea- 
sure corresponds to a probability distribution which, if evolved into the future using 
the equations of motion, remains invariant (cf. Lichtenberg and Licbcrman 1992). 
Applied to galactic dynamics, it corresponds to a time-indcpcndcnt phase space dis- 
tribution for stochastic orbits moving in some specified galactic potential. 

There are two important features about the invariant measure which are rele- 
vant in the present context. The first is that, for many time-independent potentials, 
generic ensembles of initial conditions corresponding to stochastic orbits will evolve 
on relatively short time scales towards an invariant measure. An initially localised 
ensemble of stochastic orbits will disperse and eventually, via a sort of phase mix- 
ing, evolve towards a distribution which fills all of the accessible phase space, with a 
relative weight for different phase space regions that is given by the invariant mea- 
sure. The second feature is that the invariant measure defines the natural ensemble 
of stochastic orbits to consider if one is interested in time-independent equilibrium 
configurations. In particular, it constitutes the natural building block for the con- 
struction of self-consistent models which include stochastic orbits. 

If one is interested in models of galaxies that do not exhibit strict spherical or 
axisymmetry, it is in general impossible to construct analytic equilibria. For that rea- 
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son, one is typically forced to proceed numerically via Schwarzschild's (1979) method, 
or some variant thereof. The idea underlying this method is straightforward, at least 
in principle. What one must do is (1) specify some time-independent gravitational 
potential (2) generate a large number of orbits in that potential, and then (3) 
select an ensemble of orbits which yields a mass distribution p that generates $ self- 
consistently as a solution to the Poisson equation, = AnGp. The basic point now 
is clear. If the self- consistent configuration is to constitute a true equilibrium, it must 
have the property that, when evolved into the future using the self-consistent equa- 
tions of motion, its form remains unchanged. However, when considering stochastic 
orbits, the only way to insure that this be true is to demand that the orbits be chosen 
to sample the time-independent invariant measure. 

As a concrete example, consider a nonintegrable two degree of freedom system 
characterised by a time-independent Hamiltonian H. Because H is independent of 
time, the energy E of any given orbit is conserved, so that evolution is restricted to 
a three-dimensional hypersurface of constant E. Suppose further that there exist no 
additional conserved quantities. Generically, this constant energy hypersurface will 
then contain both islands of regular orbits and a surrounding sea of stochastic orbits. 
The stochastic regions may be divided into disjoint regions by invariant tori. Suppose, 
however, that the stochastic regions are all connected, so that any given orbit can 
access all of the stochastic portions of the constant energy hypersurface. The evolution 
towards an invariant measure then corresponds to the fact that a generic ensemble of 
initial conditions, corresponding to stochastic orbits of energy E, will evolve towards 
an invariant distribution r{E), the form of which depends only on the energy and is 
independent of all other details. 

The two basic questions should now be clear: How rapid and efficient is this evo- 
lution towards the invariant measure, and how should one characterise the stochastic 
properties of this invariant measure? Answering the first of these questions involves 
probing the rate at which a generic ensemble evolves towards an "equilibrium." An- 
swering the second involves probing the properties of that equilibrium. 

At least for some model potentials, one knows that generic ensembles of initial 
conditions corresponding to stochastic orbits can evidence a very rapid coarse-grained 
evolution towards an invariant measure. For example, Kandrup & Mahon (1994) 
have shown that, for the sixth order truncation of the Toda lattice potential, generic 
ensembles of fixed energy E evolve exponentially towards an invariant measure on 
a time scale which, in physical units, is substantially shorter than tn- Moreover, 
they have shown that the rate at which ensembles of energy E evolve towards the 
invariant measure r{E) is directly related to the value of the Liapounov exponent 
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X{E). There is no guarantee that a similar behavior will be observed for all potentials. 
In particular, there may be other systems in which an invariant measure does not 
even exist. However, in such cases it would seem impossible to construct a time- 
independent self-consistent model which incorporates stochastic orbits. 

Given the importance of the invariant measure F, it is natural to search for use- 
ful characterisations of the statistical properties of orbits that sample F. One way 
in which to characterise this measure is to probe the overall degree of instability 
exhibited by orbits within the ensemble on physically relevant time scales At <tH- 
Precisely this information is contained within a calculation of short time Liapounov 
characteristic numbers 

X(At)^ Um (2) 

^ ^ 5r-(o)^o At \ 5r{0) J ^ ' 

Given the values of X{At) for all the orbits in the ensemble, one can of course 
compute A'"(X(At)), the distribution of Liapounov characteristic numbers, which de- 
pends on both the form of the invariant measure, and hence the energy, and on the 
duration of the sampling interval. One would expect that the form of this distri- 
bution should depend on the invariant measure since the avlues of the Liapounov 
exponents depend on the energy. Moreover, the distribution must depend on the 
sampling interval At since, as one samples for progressively longer times, one must 
eventually converge towards a distribution that is infinitely sharply peaked about the 
conventional Liapounov exponent. 

Even though this prescription differs substantively from the standard approach 
involving asymptotic orbital dynamics, the overall conclusions may not be all that 
different. The basic reason for this is that time averages computed for an individual 
orbit can coincide with ensemble averages generated from a collection of orbits, pro- 
vided that the ensemble is so chosen as to sample the invariant measure. This would 
imply in particular that the mean Liapounov characteristic number % associated with 
the distribution A^(X(At)) should coincide with the standard Liapounov exponent Xi 
as calculated in the usual way from a single orbit for very long times. As discussed 
more carefully in Sections 3-5, this is actually true for at least one model potential, 
namely the truncated Toda potential (cf. Kandrup & Mahon 1994). 

More generally, one might expect that the information required for the standard 
calculation of the Liapounov exponent x via a long time integration actually contains 
within it the same information as is contained within the short time distribution 
N{X{At)). Given a sequence of estimates {x(t;)}, (?' = 1,2....), computed at fixed 
intervals tj+i — ti = At, one can extract a sequence of short time Liapounov charac- 
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teristic numbers x(Aij) via the obvious prescription 

xiAti) = x{ti + At)[U + At]-xiU)U ^ 

The critical point then is that, if time and ensemble averages agree, the distribution 
of Liapounov characteristic numbers x{Ati) generated in this way should coincide (to 
within statistical errors) with the distribution N{X{At)). 

In the remaining sections of this paper, this general picture will be confirmed 
in detail for one specific model potential, namely the sixth order truncation of the 
Toda lattice. What this means is that the standard Liapounov exponent X, defined 
formally as a t — > cxo limit, actually has physical meaning on short time scales as 
well. X measures the mean stochasticity for an ensemble of orbits that samples the 
invariant measure. This implies that earher work (cf. Contopoulos & Barbanis 1989) 
tracking the time evolution of Liapounov characteristic numbers for long times can 
be exploited to extract information about the short times distribution N{X{At)). 

It is also clear that the Liapounov exponent X is only one moment of N{X{At)), 
which can also be supplemented by a consideration of higher moments. For example, 
a calculation of the dispersion ax gives a measure of the degree to which individual 
Liapounov characteristic numbers deviate from the mean. Moreover, the overall shape 
of the distribution can provide information about "stickiness" associated with islands 
of regular orbits. More precisely, one knows that stochastic orbits can sometimes get 
trapped temporarily in the neighborhood of an island, and one might expect that, 
if trapped in that neighborhood, X{At) will be substantially smaller than the mean 
value X- 

The equivalence of time and ensemble averages discussed above means that one 
can hope to apply at least one basic idea from ergodic theory to systems which 
contain both regular and stochastic orbits. Ergodic systems are systems in which 
the entire phase space is stochastic, and in which time and ensemble averages agree, 
provided that the ensemble used to construct the average is chosen appropriately 
(cf. Lichtenberg and Lieberman 1992). Thus, e.g., the standard ergodic hypothesis 
underlying equilibrium thermodynamics implies that, for an isolated Hamiltonian 
system, the natural ensemble corresponds to a microcanonical distribution, i.e., a 
uniform sampling of the constant energy hypersurface. Realistic galactic potentials 
do not seem to define ergodic systems, since even the most chaotic potentials which 
have been envisioned contain at least some regular orbits, and are characterised by 
an invariant measure which is not microcanonical (cf. Kandrup & Mahon 1994). 
However, despite deviating from ergodicity in these respects, they do appear to satisfy 
the important condition that, for the stochastic orbits, time and ensemble averages 
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agree, provided that one works with ensembles that sample the invariant measure. 

When considering regular orbits, it is physically well motivated to consider the 
properties of individual trajectories since, typically, infinitesimal changes in initial 
conditions will not lead to gross qualitative changes in the subsequent evolution. 
Thus, in particular, initially nearby trajectories remain nearby and the overall shape 
of an orbit is typically stable towards small perturbations. However, for the case 
of stochastic orbits this is no longer true. The fact that these orbits have positive 
Liapounov exponent implies that two initially nearby trajectories will typically diverge 
exponentially. In this case, it is arguably more appropriate to follow the evolution of 
an initial phase space element Az(to) than the evolution of a single material point. 
The obvious fact then is that the distribution N{X{At)) provides a characterisation 
of the average divergence of nearby trajectories and, as such, information about the 
subsequent spreading of an initially localised phase space element. 

All of this has focused on the role of Liapounov characteristic numbers in charac- 
terising the behavior of "equilibrium" ensembles that sample the invariant measure. 
However, one can also consider short time Liapounov characteristic numbers in the 
context of the approach towards equilibrium. Numerically, realisations of the invariant 
measure can be constructed by choosing any "random" ensemble of initial conditions 
corresponding to stochastic orbits and evolving these initial data forward in time until, 
statistically, they approach a samphng of a time-independent distribution. However, 
while evolving these initial conditions into the future, one can simultaneously com- 
pute short times Liapounov characteristic numbers X{At), which can provide a useful 
statistical characterisation of the overall stochasticity of the initial conditions as they 
evolve towards the invariant measure. 

3. EQUIVALENCE OF TEMPORAL AND ENSEMBLE AVERAGES 
AND THE DISTRIBUTION OF SHORT TIME X'S 

The sixth order truncation of the Toda lattice corresponds to a Hamiltonian system 
of the form 



This defines a two degree of freedom system in terms of the canonical pairs {x, Px} and 




(4) 



where 



+ '^xY - + Ix^ + xY + l^Y + ^y'- 



(5) 



8 



{y,Py}. The potential V{x,y) can be derived formally from the true Toda potential 
Vi^, y) = ^{exp[2(3)V2a; + 2y] + exp[ - 2(3)^2^; + 2y] + exp(-4y)| - ^ (6) 

by a Taylor series in x and truncated at sixth order. 

The Hamiltonian system defined by the true Toda potential (5) is integrable, but 
it is believed that all truncations at third and higher order are nonintegrable. In 
particular, one knows that the sixth order truncation is nonintegrable, and, that, 
for energies above a critical value E ^ 0.80, there exist both regular and stochastic 
orbits (cf. Contopoulos and Polymitis 1987). The Toda lattice and its truncations 
are all special in the sense that they manifest a discrete 27r/3 rotational symmetry. 
One knows (cf. Udry & Pfenniger 1988) that special symmetries of this form can 
have important effects on the form of the regular orbits and the relative abundance of 
regular and stochastic orbits. However, experience would suggest that, even though 
breaking this, or any other, symmetry will typically increase the relative abundance 
of stochastic orbits, the behavior of the stochastic orbits in this potential should 
be qualitatively similar to stochastic orbits in other potentials with different or less 
symmetry. 

Having chosen to consider a truncated Toda potential, it were perhaps appropriate 
to justify the particular truncation at sixth order. This is in fact easily done: If one 
is concerned with the evolution towards an invariant measure, one must consider a 
potential in which the orbits are confined to a compact phase space region. This is 
achieved by demanding that V{x, y) — > oo as x and y — > oo, which requires an even 
order truncation. One also wishes to consider the lowest order truncation possible, 
since this will minimise the time required in effecting the numerical computations. 
The second order truncation is of course integrable, and hence unacceptable. The 
fourth order truncation is not integrable, as proved by Yoshida et al (1988) and 
demonstrated explicitly by Udry and Martinet (1990). However, albeit nonintegrable, 
this lower order truncation exhibits relatively little chaos, perhaps because the leading 
order quartic terms in the potential (as well as the quadratic terms) are strictly 
axisymmetric. 

In order to study stochastic orbits in a given potential, one must first ascertain 
the location of the stochastic regions as a function of energy E. For sixteen different 
values of energy between E = and E = 200, surfaces of section were generated, 
plotting y and Py at successive points where randomly chosen orbits pass through 
the value x = 0. This is a useful choice of section because the Toda potential is 
symmetric under a reffection x — ^ — ,t, so that each orbit must repeatedly intersect 
the X' = hyperplane. These surfaces of section were then used to select initial data 
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corresponding to ten different stochastic orbits for each of the sixteen selected values 
of the energy. 

For each of these ten different choices of initial conditions, Liapounov exponents 
were computed in the standard way (cf. Bennetin et al 1976). This was done by 
introducing a small initial perturbation 5x — 10~^° into each orbit, continually solving 
the variational equations for the perturbation as the unperturbed orbit is evolved, and 
then renormalising the evolved perturbation back to a total amplitude 



at intervals At = 10. Successive computations of X{t) were made at the same time 
intervals, with the calculation for each orbit proceeding for a total time t = 10'^. The 
value X{t = 10'^) was then interpreted as providing an estimate for the true Liapounov 
exponent, which is of course defined only in the t — > oo limit. It was found thereby 
that, for orbits of fixed energy, the value of the Liapounov exponent is independent 
of the initial conditions, so that one can speak of a unique X{E). 

The surfaces of section were also used as an aid to construct three localised ensem- 
bles of initial conditions of fixed energy E within the stochastic phase space regions. 
Each cell was constructed as follows: One started by selecting a point {y,Py} on the 
X — surface of section which is displaced significantly from any large islands of 
regular orbits. This point was then used to define the center of a cell of specified size 
Ay and Ap^. A uniform sampling of this cell via a rectangular grid served to select 
400 pairs {y,Py}, and these pairs were then used to generate an ensemble of initial 
conditions {x,y,Px,Py}, setting x — and 



For energies E > 50, Ay — 0.2 and Apy — 0.72. For lower energies. Ay and Apy were 
chosen to be one half or one quarter as large. 

As described elsewhere (Kandrup & Mahon 1994), when these ensembles of initial 
conditions were evolved into the future, they were found to evidence a coarse-grained, 
exponential approach towards an invariant measure on a time scale t ^ 100. For this 
reason, the collections of 400 phase space coordinates at time t = 100 were taken as 
constituting random samplings of the invariant measure T{E). Short time t = 100 
Liapounov characteristic numbers were computed for these random samplings of the 
invariant measure, with X{t) for each orbit again being recorded at At = 10 intervals. 
The resulting distribution of Liapounov characteristic numbers, N{X{t = 100)), was 
then analysed to extract the first and second moments, xi^ = 100) and (^xi^ — 100). 

The first striking result derived from such an analysis is that the mean value x 
associated with the short time Liapounov characteristic numbers coincides, to within 



5z = i5x^ + 5y^ + 5pl + Splf/'' = 10 



1-10 



(7) 
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statistical uncertainties, with the long time estimate of the Liapounov exponent X. 
This is illustrated in Fig. 1, which exhibits the estimated values of the Liapounov ex- 
ponents computed in both ways. To generate this Figure, the ten long time estimates 
at each energy E were first averaged together to obtain a mean Liapounov exponent. 
The values obtained thereby were then plotted as small triangles, with error bars 
that represent the standard deviations associated with the means. The three larger 
diamonds for each value of E represent calculations of the mean values x associated 
with the short times distributions, N{X{t — 100)), generated from ensembles of 400 
orbits. 

The second significant result derived from this analysis is that the long time calcu- 
lation of X for a single stochastic orbit actually contains within it the same information 
as the distribution of short time Liapounov characteristic numbers derived for an en- 
semble of orbits. To extract this information, the values of xit) for a single orbit 
were partitioned as in Eq. (3) into intervals At = 100 to extract a different set of 
short time estimates X(Atj). The key point then is that the distribution of X(Atj)'s 
generated in this way is almost identical to N{X{At), the distribution of short time 
X's generated from the ensembles of orbits that sample the invariant measure. 

Fig. 2 illustrates the forms of these distributions for four different values of energy, 
namely E = 150, 75, 30, and 20. These distributions were generated by binning the 
data into intervals 5X — 0.05 and normalising the resulting binned distribution so 
that the most populous bin is assigned the value N{X) = 1. 

The first obvious point to be inferred from this Figure is that, as noted already, 
except for some small differences at very low values of X the two distributions are 
essentially the same. The second point is that, at least at relatively high energies, 
both distributions are extremely well approximated by a Gaussian form. However, for 
energies below a value E Ki 75 — 100, the distribution begins to acquire a statistically 
significant short X tail. When the energy is decreased further to a value as small as 
E Ki AO — 60, this tail then acquires a secondary peak. However, at very low energies, 
E 10 — 20, the two peaks eventually merge into a single distribution which deviates 
significantly from a Gaussian. 

The observed structure at small values of x seems to be associated with the increas- 
ing predominance of regular orbits at lower energies. It is well known (cf. MacKay 
et al 1984a,b, Lau et al 1991) that stochastic orbits that stray too close to islands 
of regularity tend to get trapped for relatively long periods of time around these is- 
lands, and it would appear that, while in these regions, their Liapounov characteristic 
numbers are substantially reduced. 

The differences between the two distribution at small values ofx < 0.1 — 0.15 arise 
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because the alleged sampling of the invariant measure is contaminated by a small 
number of regular orbits. In constructing these samplings, one selected an initial 
phase space region which seemed, at least superficially, to be comprised completely 
of stochastic orbits. However, these stochastic regions actually contain a nonzero 
measure of regular orbits embedded in the stochastic sea. 

The basic conclusion, therefore, is that, because of the equivalence of time and 
ensemble averages, an evaluation of X(t) for a single orbit over long times contains the 
same information as the distribution of short time Liapounov characteristic numbers 
associated with an ensemble of orbits that samples the invariant measure. 

4. LIAPOUNOV CHARACTERISTIC NUMBERS AND THE 
APPROACH TOWARDS AN INVARIANT MEASURE 

The preceding section focused on the form of the distribution of Liapounov character- 
istic numbers associated with an initial ensemble that samples the invariant measure. 
However, as discussed in Section 2, one can also consider the distribution N{X) as- 
sociated with other ensembles as well. In particular, it is natural to compute the 
distribution of Liapounov characteristic numbers for the initially localised ensembles 
considered in Section 3 as they evolve towards the invariant measure. Doing this 
enables one to understand how the evolution towards an invariant measure correlates 
with the overall stochasticity of the orbits, a problem already considered in a slightly 
different context by Kandrup & Mahon (1994). 

The results for one particular energy, namely E = 50, are illustrated in Fig. 3. 
The data plotted in this Figure were generated by selecting the three initially localised 
ensembles discussed in Section 3, and computing X{t) for each of the 400 orbits in 
each ensemble at t = 10 intervals for a total time oi t = 200. The resulting X's were 
then partitioned into collections of short time estimates X(Atj) at At = 10 intervals. 
Fig. 3 exhibits the mean x(Atj) and the associated dispersions ax{Ati) for these 
three different ensembles. 

In this Figure, the first interval At — 100 is interpreted as corresponding to the 
period during which the ensemble evolves towards the invariant measure. The phase 
space coordinates ait — 100 correspond to the sampling of the invariant measure used 
for the computations in Section 3, and the remaining At — 100 are thus interpreted 
as corresponding to the evolution of an initial ensemble that samples the invariant 
measure. 

From this Figure several conclusions are apparent. The first is that even at very 
early times, when the ensemble deviates significantly from a sampling of the invariant 
measure, the values of xi^U) and a^iAti) do not deviate all that much from the 
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"equilibrium" values associated with the invariant measure. At early times, x tends 
to be somewhat larger than the value associated with the invariant measure. This is 
easily understood by observing that the initial ensemble was located in a phase space 
region particularly far from any large islands of regular orbits where, overall, one 
might expect a higher degree of stochasticity. At early times, one also observes that 
(j^ is somewhat smaller than the value associated with the invariant measure. This 
is again easily understood by observing that the initial ensemble is constructed from 
a collection of phase space points that are relatively close together and significantly 
displaced from any large islands. 

It is also apparent that, by a time t = 100, x and ax have asymptoted towards 
constant values, this corroborating the expectation that the phase space coordinates 
at t = 100 can be interpreted at least appropriately as constituting a random real- 
isation of the invariant measure. For higher values of energy, the approach towards 
the invariant measure is even more rapid. Indeed, as discussed more extensively 
in Kandrup & Mahon (1994), there exists a direct one-to-one correlation between 
the exponential approach towards an invariant measure and the value of the Lia- 
pounov exponent. Specifically, as was demonstrated in especial detail for energies 
10 < E < 75, increasing E increases both (1) the value of the Liapounov exponent 
X{E) associated with stochastic orbits of that energy and (2) the exponential rate 
A{E) associated with the approach towards the invariant measure, in such a fashion 
that the ratio TZ{E) = A/X is approximately constant, independent of energy. The 
fact that the mean x{t) associated with the evolution of the initially localised ensem- 
ble is rather close in value to the Liapounov exponent X, as defined in a i — > oo limit, 
explains why there can exist a direct connection between the "equilibrium" X and the 
"nonequilibrium" evolution towards an invariant measure. 

Finally, it should perhaps be noted that the first points in Fig. 3 at i = 10 
may be somewhat suspect. Specifically, the values of X{t — 10) may be significantly 
infiuenced by the particular choice of initial perturbation, namely Sx — 10~^° and 
5y = 5px = 5py = 0. 

5. THE EFFECTS OF A VARYING SAMPLING TIME 

In this Section, attention focuses on how the form of the distribution of short time 
X's extracted from a single long time integration, or generated from an ensemble that 
samples the invariant measure, varies as a function of the sampling interval At. One 
anticipates physically that, as the length of the sampling interval increases, the mean 
Liapounov characteristic number x will provide an increasingly better characterisation 
of the overall stochasticity, and that the dispersion ax associated with the distribution 
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will tend to zero for — > oo. 

Given a long time calculation of X{t) for a single orbit, one can extract short time 
X(Aij)'s for various choices of sampling interval and compare the results. Fig. 4a 
exhibits the values of X(Aij) for one long time integration of an orbit with E — 150, 
for a sampling interval At — 10. It is clear that X(Aij) shows rapid large amphtude 
variability, so that the degree of stochasticity exhibited by the orbit at time t may be 
substantially different at times t ± At. However, if one considers a longer sampling 
interval, the variability is substantially reduced. This is illustrated in Fig. 4b, which 
analyses the same data for sampling intervals At — 100. This Figure was generated 
from Fig. 4a by taking time averages of successive collections of ten intervals with 
At = 10. 

The differences between Figs. 4a and 4b can be quantified by specifying the 
maximum and minimum values of X(Atj), as well as the standard deviation about 
the mean. For the shorter time sampling, the minimum and maximum values of X are 
—0.125 and 1.670, and the standard deviation is 0.2997. For the longer time sampling, 
the minimum and maximum values are 0.078 and 1.053, and the standard deviation 
is approximately half as large, namely 0.1449. However, even though the longer time 
sampling is smoother, nontrivial structures are still observed. In particular, there is 
a significant dip in the value of X between t = 4000 and 4300 which is observed in 
both figures, where X{Ati) decreases appreciably to a value < 0.1. 

To further elucidate these sorts of differences, one can also examine the form of 
the distribution N{X{At)) as a function of the sampling interval. Figures 5a and 5b 
illustrate the form of N{X{At)) for two different energies, E — 150 and 50, and three 
different sampling intervals. At — 10, 40, and 100. These figures were constructed 
by partitioning each of the ten t — 10"^ calculations into a collection of intervals of 
fixed length At, combining the data for all ten orbits, and constructing a binned 
distribution with AX = 0.05. For both values of the energy, it is clear that the 
width of the distribution decreases as the sampling interval becomes longer. It is also 
clear that, for each choice of sampling interval, the relative width of the low energy 
distribution is shghtly larger than that of the high energy distribution. For E — 150, 
the distribution is essentially Gaussian for the two larger sampling intervals, whereas, 
for £■ = 50 one sees indications of a low X tail. 

It is natural to ask whether the systematic decrease in the width of the distribu- 
tion, as probed by the dispersion, is exponential in time, or whether it is better fit 
by a power law. The answer is that, for all values of the energy, the dispersion is well 
fit by a power law a-)^ oc At~^ , where p is a positive constant. The goodness of fit is 
illustrated in Fig. 6 for three different energies, E — 50, 125, and 175. 
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One might also ask how the value p depends on the energy E. It is difficult to 
extract a precise estimate of p for any given energy because the best fit value depends 
sensitively on whether or not there exists one or two initial conditions where the orbit 
spends a large amount of time in a region of especially low X. However, one can still 
extract a best fit value of p for each energy. The result of such an analysis is exhibited 
in Fig. 7. Because of the large scatter in this plot, it is difficult to extract a precise 
functional form for p{E). However, two trends are unambiguous: First, it is clear 
that, overvall, the value of p increases with increasing energy, and second, it would 
appear that the value of p approaches 0.5 at large energy. 

A value of p = 0.5 is relatively easy to explain. The distribution A^(X(At)) for 
large At can be viewed as a convolution of a large number of distributions A^(X(Atj)) 
for shorter time intervals At = 10. Suppose however, that the values of X(Atj) for 
successive At = 10 intervals are completely uncorrelated. In this case, the central 
limits theorem guarantees that the long time distribution will converge towards a 
Gaussian with a dispersion that scales as k^^^"^, where k denotes the total number 
of short time intervals. The fact that this number grows linearly in time would then 
imply a dispersion that scales as t~^/^. If these successive intervals are not completely 
uncorrelated, then the dispersion should decrease more slowly. Fig. 7 therefore 
indicates that, at low energies, there is substantial correlation between successsivc 
At = 10 intervals, but that this correlation decreases for higher energies. 
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Fig. 1. Liapounov exponents as a function of energy E estimated in two different 
ways. The small triangles represent x, calculated in the usual way by following ten 
different orbits for a total time t — 10^. The error bars represent the associated 
standard deviations. The larger diamonds represent calculations of the mean x{t) for 
ensembles of 400 orbits over a period t — 100. 

Fig. 2. (a) The binned distribution of Liapounov characteristic numbers x(At) 
for At — 100, calculated in two different ways for orbits with E — 100. The solid 
curve represents the distribution obtained from an initial ensemble that sampled 
the invariant measure. The dashed curve represents the distribution obtained from 
partitioning the curve x{t) for t — 10000 into 100 segments of length At — 100, and 
construction x(Atj) for each segment, (b) The same for E — 75. (c) The same for 
E = 30. (d) The same for E = 20. 

Fig. 3. The mean x(Ati), and the associated dispersion ax{Ati), computed at 
At — 10 intervals for a total time t — 200 for three initially localised ensembles of 400 
stochastic orbits with E — 50. The evolved phase space coordinates t — 100 were 
interpreted as constituting random realisations of the invariant measure, and used as 
initial data for the calculations described in Section 3. The solid lines represent the 
values of x and associated with the invariant measure. 

Fig. 4. The Liapounov characteristic number x(^) for one orbit with E = 100 
partitioned into a collection of short times estimates x(At), with (a) At = 10 and (b) 
At = 100. 

Fig. 5. (a) The distribution of short times x(Ai) for E — 150, extracted from 
the long time xit)- The solid curve represents a sampling interval At = 100, the 
dot-dashed line At — 40, and the dashed hne At — 10. (b) The same for E — 50. 

Fig. 6. The dispersion as a function of samphng interval At for energies £■ = 50 
(bottom curve), 125 (middle curve), and 175 (top curve). The sohd curves represent 
least squares power law fits ax{At) oc At~^. The curves for £■ = 50 and E — 175 
were displaced respectively upwards and downwards by crx = 0.125 to provide a less 
cluttered diagram. 

Fig. 7. The exponent p associated with the least squares fit ax oc At~^', plotted as 
a function of energy E. 
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